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Abstract 

The  numerical  simulation  of  steady  planar  two-dimensional,  laminar  flow  of  an 
incompressible  fluid  through  an  abruptly  contracting  channel  using  spectral  domain 
decomposition  methods  is  described.  The  key  features  of  the  method  are  the 
decomposition  of  the  flow  region  into  a  number  of  rectangular  subregions  and 
spectral  approximations  which  are  pointwise  Cl  wontinuous  across  subregion 
interfaces.  Spectral  approximations  to  the  solution  are  obtained  for  Reynolds 
numbers  in  the  range  10,  500).  The  size  of  the  salient  corner  vortex  decreases  as  the 
Reynolds  number  increases  from  0  to  around  45.  As  the  Reynolds  number  is 
increased  further  the  vortex  grows  slowly.  A  vortex  is  detected  downstream  of  the 
contraction  at  a  Reynolds  number  of  around  175  that  continues  to  grow  as  the 
Reynolds  number  is  increased  further. 

Research  was  supported  in  part  for  the  second  author  by  the  National 
Aeronautics  and  Space  Administration  under  NASA  Contract  No.  NAS1-18605  while  he 
was  in  residence  at  the  Institute  for  Computer  Applications  in  Science  and 
Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23665. 


L  introduction 

The  numerical  simulation  of  the  flow  in  a  constricted  channel  is  considered 
using  spectral  domain  decomposition  techniques.  The  first  numerical  solutions  to 
this  problem  were  obtained  by  Dennis  and  Smith  (1980)  using  a  finite  difference 
discretization  of  the  stream  function  -  vorticity  formulation  of  the  governing 
equations.  They  did  not  detect  a  downstream  recirculation  region  caused  by  the 
flow  separating  at  the  corner  despite  using  a  very  fine  uniform  grid.  In  a  recent 
paper  Hunt  (1989)  uses  a  non-uniform  grid  to  ensure  a  dense  distribution  of  grid 
points  in  the  regions  of  the  flow  which  need  to  be  resolved.  He  places  a  locally 
fine  mesh  downstream  of  the  constriction  and  around  the  re-entrant  corner  since 
this  gives  rise  to  a  singularity  in  the  vorticity. 

In  this  paper  the  governing  equations  are  written  in  terms  of  the  stream 
function.  This  means  that  mass  is  conserved  identically.  The  flow  region  is 
divided  into  a  number  of  rectangular  subdomains.  In  each  of  these  subdomains  the 
stream  function  is  approximated  by  a  truncated  double  Chebyshev  expansion.  The 
expansion  coefficients  are  determined  by  collocating  the  governing  equation  and 
appropriate  interface  continuity  conditions  between  subdomains.  The  decomposition 
is  such  that  the  Chebyshev  collocation  points  are  distributed  densely  around  the 
corner  and  near  solid  boundaries.  This  is  important  computationally  in  order  to 
resolve  the  main  features  of  the  flow  efficiently. 

In  previous  work,  Karageorghis  and  Phillips  (1989a)  consider  nonconforming 
subdomains  because  of  their  ease  of  implementation.  Although  this  strategy  works 
well  for  the  Stokes  problem,  a  lack  of  interface  continuity  appears  for  the  Navier- 
Stokes  problem  (Karageorghis  and  Phillips  (1989b))  for  values  of  the  Reynolds 
number  around  200,  eventually  causing  the  method  to  break  down  completely.  The 
spectral  approximations  obtained  are  not  pointwise  continuous  across  the  subdomain 
interface.  In  this  paper  a  collocation  strategy  is  used  which  generates  pointwise  C‘ 
continuous  approximations  across  the  interfaces.  As  a  result  no  computational  limit 
on  the  value  of  the  Reynolds  number  has  yet  been  encountered.  Numerical 
solutions  are  obtained  for  values  of  the  Reynolds  number  up  to  500. 

The  matrices  arising  from  spectral  discretizations  are  not  sparse  like  their 
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finite  difference  and  finite  element  counterparts.  However,  when  used  in 
conjunction  with  domain  decomposition  techniques  these  matrices  possess  a  block 
tridiagonal  structure  which  can  be  exploited  when  designing  methods  of  solution. 
Here  we  use  a  subroutine  from  the  NAG  Library  which  solves  almost  block  diagonal 
systems  (Brankin  and  Gladwell  (1990)).  This  method  is  as  efficient  as  the 
capacitance  matrix  method  in  terms  of  cost  and  storage  but  is  found  to  be  more 
stable  (Karageorghis  and  Phillips  (1990)). 

Spectral  approximations  to  the  solution  of  this  problem  are  obtained  for 
Reynolds  numbers  in  the  range  (0,  500],  A  comparison  with  the  work  of  Dennis  and 
Smith  (1980)  shows  good  agreement  between  the  sets  of  results  in  terms  of  stream 
function  contours  in  the  bulk  of  the  flow  and  the  description  of  the  salient  corner 
vortex.  Although  they  did  not  find  a  downstream  recirculation  region  there  is  a 
hint  of  its  existence  on  their  finest  glides  for  Reynolds  numbers  above  1000.  They 
probably  failed  to  detect  this  behaviour  in  their  numerical  simulations  because  the 
grid  employed  was  not  fine  enough  in  this  region.  In  our  calculations  separation 
first  appears  at  around  Re  =  175  which  is  slightly  earlier  than  Hunt  (1989)  predicts. 
The  separation  length  and  the  strength  of  the  vortex  increase  as  Re  increases,  the 
length  being  roughly  proportional  to  Re.  The  spectral  collocation  method  predicts 
longer  separation  lengths  than  Hunt  (1989).  On  the  finest  grids  used  the  spectral 
collocation  method  employs  about  a  third  of  the  number  of  degrees  of  freedom  as 
Hunt  (1989)  and  a  fifteenth  of  the  number  of  degrees  of  freedom  as  Dennis  and 
Smith  (1980). 

2.  The  Governing  Equations 

We  consider  the  steady  laminar  flow  of  an  incompressible  fluid  through  an 
abruptly  contracting  channel  with  walls  aty-±lforx<0,  y  =  ±*  for  x  >  0 

A* 

and  1/2  <;  |yl  <(  1  for  x  =  0.  Upstream  we  impose  parabolic  Poiseuille  flow  and 
we  suppose  that  the  flow  is  parabolic  again  far  enough  downstream.  Since  the  flow 
is  symmetric  about  y  =  0  it  is  only  necessary  to  seek  a  solution  for  y  ;>  0. 

In  terms  of  non-dimensional  variables  the  incompressible  Navier-Stokes 
equations  are 


(v  .  7)v  =  -  7p  +  (Re)'1  72v 


(2.1) 
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7  .  v  =  0 


12.2) 


where  v  =  (u,  v)  is  the  velocity  vector,  p  is  the  pressure  and  Re  is  the  Reynolds 
number.  The  introduction  of  a  stream  function,  0(x,  y),  defined  by 


u 


30  _  30 

3y  V  —  3x 


(2.3) 


I 

4 


means  that  the  continuity  equation  (2.2)  is  satisfied  identically  and  (2.1)  becomes, 
after  elimination  of  the  pressure. 


"0  -  Re 


~  |-  (7-0)  -  |-(7'0) 

3 y  3x  3x  3y 


=  0 


The  boundary  conditions  are 


(2.4) 


0  =  1,  on  y  =  1,  x^O  and  y  =  ^  ,  x  ^  0,  (2.5) 

0  =  1,  ||  =  0  on  x  =  0,  \  £  y  <;  1  ,  (2.6) 

0=0,  =  0  on  y  =  0.  (2.7) 

3y 


The  conditions  on  y  =  0  preserve  the  symmetry  of  the  flow  while  the  other 
conditions  represent  the  no-slip  conditions. 

The  governing  equation  (2.4)  is  a  nonlinear  partial  differential  equation  for 
the  stream  function.  It  is  solved  in  an  iterative  fashion  using  a  Newton-type 
method  to  linearize  it  (Phillips  (1984)).  It  is  the  linearized  partial  differential 
equation  which  is  solved  at  each  Newton  step  using  a  spectral  collocation  method. 
Let  us  write  (2.4)  in  operator  form, 

1(0)  =  0  (2.8) 
where  L  is  a  nonlinear  operator.  Suppose  that  0*  is  some  approximation  to  the 
solution  of  (2.8).  To  obtain  an  improved  approximation  wr  replace  L  by  its 
linearization  about  0*  and  then  solve  the  problem 

L'(0*)d  =  -  L(0*)  ,  (2.9) 
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where  L'($)  is  the  Frechet  derivative  of  L  at  t£>  defined  by 


L'(m  =  ”*<t>  -  Re 


~  (V-0)  - 

3y  3x 


3x  3y 


V7-H»  ~ 
3x  3y 


3y  3x 


(2.10) 

The  new  approximation  to  the  solution  is  thus  H>*  -r  <t> •  This  completes  a  single 
Newton  step  and  is  repeated  until  convergence  is  reached. 

3.  Domain  Decomposition  and  Spectral  Approximation 

Since  spectral  methods  are  most  easily  applied  to  problems  defined  in 
rectangular  regions,  the  natural  domain  decomposition  of  the  flow  region  comprises 
three  rectangular  subdomains  with  common  point  (0,  1/2)  as  shown  in  Figure  1.  Such 
a  domain  decomposition  may  be  termed  conforming  since  the  sides  of  each  rectangle 
are  contiguous  with  either  a  boundary  segment  or  a  side  of  another  rectangle.  A 
further  advantage  of  this  decomposition  is  that  it  is  possible  to  construct  a 
collocation  scheme  which  results  in  pointwise  C1  continuous  approximations  across 
the  subdomain  interfaces. 

The  flow  region  is  truncated  upstream  and  downstream  of  the  constriction  at 
finite  distances  ht  and  h;  from  the  origin,  respectively.  The  distances  h,  and  h; 
need  to  be  sufficiently  large  so  that  the  flow  is  fully  developed  in  the  entry  and 
exit  sections.  The  domain  truncation  means  that  additional  boundary  conditions 
need  to  be  imposed  on  entry  and  exit,  namely  that  the  normal  derivative  of 
vanishes  on  entry  and  exit. 


B  C(0,1) 


1 

D(0,  1/2) 

II 

III 

A(-hj,  0)  GCO,  0)  F(h2,  0) 

Fig.  1  Three  subdomain  decomposition  of  flow  region. 
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In  each  subdomain  the  stream  function  tf>(x,  y)  is  approximated  by  ?0k(x,  y), 

where 

Mk  Nk 

tfk(x,  y)  =  Gk(y)  +  ^2  XI  Wm(x)  P"(y)  ’  k  =  l’  H’  ,n*  {3A) 

m-Mg  n-Nk 

and 

G!(y)  -  Gn(y)  =  iy(3-y-)  ,  GU1(y)  =  Gl(2y)-, 

4m 

M0  “  N0  -  M0  -  No  -  Mo“  -  No"  -  4  • 

k  k 

The  polynomials  (Wm(x)}  and  (Pn(y)}  are  modified  shifted  Chebyshev  polynomials 
which  satisfy  identically  all  the  boundary  conditions  with  the  exception  of  the 
conditions  along  the  vertical  wall  CD.  For  example,  in  region  I  we  have 

Wjnlx)  =  Tjn(x)  +  alm  t[(x)  +  ^T^x)  ,  2  £  m  £  M1  ,  (3.2) 

where  t{„(x)  .  0  <;  m  ^  ,  are  shifted  Chebyshev  polynomials  on  [— h^,  0] 

defined  by 

r  r2x  -(-  hu 

T}n(x)  =  Tm(  — — ^)  , 

and  a{„  and  /?{„  are  given  by 

«m  =  (-l)m  nr  ,  0{m  =  (-l)m(nr  -  1)  ,  2  £  m  £  M1  . 

Similarly,  we  can  show  that 

pH  (y)  =  f  Uy)  +  d*T  J(y)  +  fiKv [0(y)  ,  2  £  n  <;  N1  ,  (3.3) 

where  T J  (y)  ,  0  ^  n  <;  ,  are  the  shifted  polynomials  on  (1/2,  1]  defined  by 


T‘(y) 


Tn  (4y-3) 


and  ai  ,  q\  are  given  by 

a1,  =  —  n:  —  1,  =  n;  —  1  . 

The  basis  functions  in  regions  II  and  III  are  defined  in  a  similar  fashion. 


4.  Collocation  Strategy 

The  coefficients  (a^n)  ,  k  =  I,  II,  III  must  be  determined  for  the  spectral 
approximations  (3.1)  to  be  defined  everywhere  in  the  truncated  flow  region. 
Suppose  that  at  the  beginning  of  a  Newton  step  we  have  approximations  a*„  to  the 
expansion  coefficients  of  iC^(x,y).  The  numerical  solution  of  the  linearized  equation 
(2.9)  determines  the  expansion  coefficients  (tfa)™  of  0k,  the  correction  to  lbk,  in 
region  k  for  k  =  I,  II,  III.  We  describe  in  some  detail  the  process  of  calculating  the 
coefficients  (5a)™  Once  these  are  known,  the  updated  approximation  to  the  stream 
function  is  just  #*k  -+-  0k  ,  which  can  be  found  by  simply  adding  together  a*k  and 

(6a)™  for  each  value  of  m,  n  and  k. 

The  linearized  partial  differential  equation  (2.9)  is  collocated  in  each 
subdomain  and  the  approximations  in  the  three  subregions  are  patched  by  imposing 
the  correct  order  of  weak  continuity  across  the  subregion  interfaces.  The 
collocation  points  in  a  subdomain  are  chosen  to  be  those  points  at  which  the 
Chebyshev  polynomial  of  highest  degree  used  in  the  representation  in  that 
subdomain  attains  its  extreme  values.  This  choice  gives  rise  to  optimal 

approximation  properties  of  smooth  functions.  It  can  also  be  shown  that  when 
Gaussian  quadrature  rules  based  on  these  points  are  used  to  evaluate  the  integrals 
appearing  in  Galerkin  formulations  of  certain  differential  equations,  then  the 
resulting  equations  are  equivalent  to  those  determined  by  collocating  the 
differential  equation  at  the  same  set  of  points.  Thus,  in  region  I,  for  example,  the 
collocation  points  are  given  by 


hj(Xj  —  1) 
2 


yj  +  3 


4 


(4.1) 


where 
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cos  (i-^)  ,  0  <;  ,  <  M1 

M1 

cos  (^)  ,  0  £  j  <;  N1  . 

N1 


Boundary  conditions 

Due  to  the  choice  of  modified  Chebyshev  polynomials  as  trial  functions  in 
the  expansions  (3.1)  the  boundary  conditions  are  automatically  satisfied  except  along 
x  =  0  (1/2  ^  y  <;  1)  in  region  I.  Along  this  portion  of  the  boundary  there  are 
N'^  +  l  collocation  points.  We  deduce  from  (3.1)  that  #  and  are  polynomials  of 
degree  in  y  along  x  =  0  (1/2  ^  y  £  1)  each  depending  on  —  1  degrees  of 

freedom.  Therefore  collocation  of  the  boundary  conditions  at  the  N*  —  1  distinct 
points  (0,  y^j),  j  =  0,  1,  ...,  —  2,  ensures  that  the  no-slip  boundary  conditions 

are  satisfied  identically  along  this  part  of  the  boundary. 

Interface  continuity  conditions 

We  impose  C3  continuity  of  the  stream  f  unction  across  the  interfaces 
between  the  subdomains.  These  conditions  are  collocated  in  such  a  way  so  as  to 
yield  approximations  that  are  pointwise  Cl  continuous  across  the  interfaces,  but 
whose  second  and  third  normal  derivatives  are  continuous  only  at  the  interface 
collocation  points.  If  the  initial  approximation  to  the  solution  of  the  problem 
satisfies  these  conditions  then  these  continuity  conditions  are  imposed  on  <t>  at  each 
Newton  step. 

Let  us  examine  the  interface  y  =  l/2(  —  h,  <,  x  <,  0)  between  subregions  I 
and  11,  for  example.  Along  y  =  1/2,  ^  and  <t>^  are  polynomials  of  degree  M  each 
possessing  —  1  degrees  of  freedom.  We  collocate  at  a  sufficient  number  of 
points  on  the  interface  to  ensure  that  ^  —  <t> ^  is  identically  zero.  This  is 
achieved  if  we  impose 

0I(xJ  ,  i)  -  4>U(x\  ,  i)  =  0  ,  i  -  2,  3,  ...,  M1  -  2,  (4.2) 

and,  in  addition, 

*II(0,  =  0>  IT  (°’  2}  = 


0 


(4.3) 
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Similarly,  continuity  of  <t> y  is  obtained  by  imposing 


§*  (xl  l)  _  ^(xi  L) 

3y  1  »  ’  2J  3y  '  i  ’  2 


=  0  ,  i  =  2,  3 . M1 


(4.4) 


and 


3  <t>u 
3y 


(0,  b  =  0 


3x3y 


=  0 


(4.5) 


The  conditions  which  represent  the  continuity  of  the  second  and  third  normal 

derivatives  of  0  are  collocated  at  the  —  3  points  (x{  ,  1/2)  ,  i  =  2,  3 . 

M*—  2.  Thus,  these  derivatives  are  not  pointwise  continuous  across  the  interface. 
Moffatt  (1964)  shows  that  the  leading  term  in  the  singular  expansion  of  the  Stokes 
solution  about  the  re-entrant  corner  (0,  1/2)  is  O(r^)  where  X  1.5445  and  so  it 
would  be  inconsistent  to  impose  pointwise  continuity  of  these  higher  derivatives 
across  the  interface. 

The  same  collocation  strategy  is  applied  across  the  interface  x  =  0 
(0  ^  y  <1  1/2)  between  subregions  II  and  III.  As  a  result  the  correction  <t>  and  its 
normal  derivative  are  pointwise  continuous  across  the  interfaces.  The  updated 
approximations  to  the  stream  function  also  possess  the  same  order  of  pointwise 
continuity. 


Linearized  differential  equation 

The  linearized  differential  equation  (2.9)  is  collocated  in  each  subregion  at  all 
points  on  the  collocation  grid  with  the  exception  of  those  on  or  one  in  from  each 
subregion  boundary,  i.e., 

(xk  ,  yk)  ,  i=2,  ....  Mk  -  2,  j  =  2 . Nk  -  2, 


for  k  -  I,  II  or  III. 


When  the  spectral  collocation  equations  resulting  from  the  boundary 
conditions,  interface  continuity  conditions  and  differential  equation  are  added 
together,  they  yield  a  total  of  [(M^  —  1)  (N^  —  1)  4-  (M^  —  1)  (N^  —  1)  -f 
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(M^  —  1)(N^  —  3)1  linear  equations  which  is  equal  to  the  number  of  unknown 
expansion  coefficients  (5a))™  for  k  =  I,  II,  III.  Therefore,  provided  the  coefficient 
matrix  is  non-singular,  this  system  of  equations  possesses  a  unique  solution. 


5^  Direct  Method  of  solution 

The  global  system  for  the  expansion  coefficients  is 

Gz  =  r  (5.1) 


where 


A 

0 

0 

(5a)1 

B 

0 

0 

c 

0 

.  Z.  = 

(5a)11 

0 

D 

0 

0 

E 

_ 

. 

l3 


r-4 


-5 


(5.2) 


The  vectors  ($a)^,  (<5a)^  and  (6a)^  contain  the  unknown  coefficients  in  the 
expansions  in  subregions  I,  II  and  III,  respectively.  The  first,  third  and  fifth  blocks 
of  rows  in  (5.2)  correspond  to  the  collocation  of  the  linearized  differential  equation 
and  boundary  conditions  not  satisfied  by  the  approximations  in  subregions  1,  II,  and 
III,  respectively.  The  second  and  fourth  blocks  of  rows  correspond  to  the 
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collocation  of  the  interface  continuity  conditions  between  subregions  1  and  11,  and  II 
and  III,  respectively. 

The  system  (5.1)  is  solved  using  an  application  of  a  standard  production  code 
designed  for  the  solution  of  almost  block  diagonal  systems  (Brankin  and  Gladwell 
(1990)).  In  a  recent  survey  of  direct  methods  for  solving  systems  of  equations 
arising  from  spectral  domain  decomposition  methods,  Karageorghis  and  Phillips  1 1 990) 
found  this  solver  to  be  superior  to  other  techniques  with  regard  to  cost,  stability 
and  storage.  The  code  uses  a  modified  column  elimination  procedure  with  alternate 
row  and  column  pivoting  based  on  an  algorithm  originally  described  in  Varah  (1976) 
and  Diaz,  Fairweather  and  Keast  (1983)  and  is  intended  to  solve  systems  of  the  fonm 
shown  in  Figure  2.  These  systems  comprise  rectangular  blocks  along  the  diagonal 
and  are  such  that  no  three  successive  blocks  have  columns  in  common. 


Fig. 2  Almost  block  diagonal  form. 

The  global  spectral  collocation  matrix  G  in  (5.1)  may  be  written  in  almost 
block  diagonal  form  in  the  obvious  way: 
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(5.3) 


The  matrix  G  has  five  non-zero  blocks,  namely, 


p2  Q: 

R  *5  S-> 

N  ■ 

Qi  V 

•  N  ■ 

Sl  T,  j 

and  jT->j 


Unfortunately  the  form  (5.3)  is  not  of  the  structure  required  by  the  almost  block 
diagonal  solver  due  to  too  much  overlap  between  blocks  2  and  3,  and  3  and  4. 
However,  the  transpose  of  G  is  of  the  required  form,  with  three  blocks,  namely, 


qT 


,T,  .T 


and 


T 

One  may,  therefore,  decompose  the  transpose  of  the  global  matrix  G  using  the 

existing  NAG  subroutine  F01LHF  and  subsequently  solve  using  the  transpose  of  the 

T  -  t 

decomposed  form  of  G  ,  say  G  ,  the  system 


(5.4) 


with  the  NAG  subroutine  F04LHF.  Further  details  of  the  implementation  may  be 
found  in  Karageorghis  and  Phillips  (1990). 
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6.  Numerical  Results 

The  results  of  numerical  calculations  of  laminar  flow  through  a  3:1 
constricted  channel  are  presented  and  comparisons  made  with  the  work  of  Dennis  and 
Smith  (1980)  and  Hunt  (1989).  Dennis  and  Smith  (1980)  use  a  finite  difference 
discretization  of  the  stream  function-vorticity  formulation  of  the  Navier-Stokes 
equations  on  a  uniform  grid.  An  upwind  differencing  scheme  of  Dennis  and  Hudson 
(1978)  is  used  to  approximate  the  vorticity  transport  equation.  This  is  essential  for 
reasonably  high  values  of  Re  in  order  to  maintain  the  diagonal  dominance  in  the 
approximating  sets  of  difference  equations.  The  loss  of  diagonal  dominance  presents 
difficulties  in  the  convergence  of  iterative  solution  techniques.  If  upwmding  is  not 
used  then  diagonal  dominance  can  only  be  established  either  by  sufficiently 
decreasing  the  mesh  size,  which  may  be  prohibitively  expensive,  or  by  using  a  non- 
uniform  grid.  Hunt  (1989)  also  uses  a  finite  difference  discretization  of  the  stream 
function-vorticity  formulation  but  on  a  non-uniform  grid.  The  vorticity  unknowns 
are  eliminated  to  give  a  system  solely  in  terms  of  unknown  values  of  the  stream 
function.  Hunt  (1989)  investigates  the  application  of  artificial  viscosity  or 
upwmding  but  contrary  to  expectations  finds  that  the  scheme  with  upwinding  fails 
to  converge  for  Re  >  500  even  though  one  would  expect  the  addition  of  a  viscous’ 
like  term  to  have  a  stabilizing  effect.  Further  for  a  given  value  of  Re  convergence 
becomes  increasingly  more  difficult  on  successively  finer  grids.  The  opposing 
conclusions  on  the  application  of  upwind  differencing  must  be  due  to  the  different 
nature  of  the  grids  used  in  these  studies.  The  ability  of  these  two  methods  and 
the  spectral  algorithm  to  describe  the  main  features  of  the  flow  is  examined. 

The  solution  of  the  Stokes  problem  (Re  =  0)  is  used  as  the  initial 

approximation  to  the  solution  of  the  Navier-Stokes  problem  for  0  <  Re  150. 

However,  for  Re  <1  150  the  Newton  process  is  robust  enough  to  converge  from  an 

initial  approximation  that  is  zero  everywhere.  For  Re  >  150  continuation  in  Re  is 

required  for  convergence  in  increments  of  Re  of  50.  This  means,  for  example,  that 
the  converged  numerical  solution  for  Re  =  150  is  used  as  the  initial  approximation 
when  calculating  the  numerical  solution  for  Re  =  200. 

The  Newton  process  is  terminated  when  the  maximum  residual  is  less  than 
IQ'7  in  magnitude.  This  invariably  occurs  after  six  Newton  iterations.  The 
expansion  coefficients  are  also  checked  for  convergence.  The  last  two  iterations 
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are  found  not  to  affect  the  first  eight  significant  digits  of  ail  the  expansion 
coefficients. 

Numerical  experiments  are  performed  on  a  number  of  collocation  grids.  This 
is  necessary  in  order  to  ensure  that  the  approximations  obtained  are  not  strongly 
grid  dependent.  The  following  discretization  parameters  define  the  collocation 
grids  used  in  this  work  and  the  total  number  of  degrees  of  freedom: 

(1)  M1  =  M11  =  24,  N1  =  20,  N11  =  N111  =  16,  M111  =  16  (977  degrees  of 
freedom); 

(2)  M1  =  M11  =  24,  N1  =  N'11  =  N111  =  20,  M111  =  24  (1265  degrees  of 
freedom); 

(3)  M1  =  M11  =  24,  N1  =  N11  =  N111  =  20,  M111  =  32  (1401  degrees  of 
freedom); 

(4)  M1  =  M11  =  24,  N1  =  N11  =  N111  =  20,  M111  =  40  (1537  degrees  of 
f  reedom). 

For  Re  <;  100  good  agreement  is  obtained  on  all  four  grids.  For  100  <;Re  <^300 
good  agreement  is  reached  on  grids  (2)  —  (4).  Thereafter,  for  Re  ^  300,  grids  (3) 
and  (4)  produced  almost  identical  results. 

The  upstream  and  downstream  truncation  lengths  are  chosen  to  be  1.0  and  3.5, 
respectively.  It  is  necessary  to  experiment  with  the  number  of  degrees  of  freedom 
in  subregion  III  since  for  an  insufficient  number  the  flow  is  not  adequately  resolved 
particularly  near  the  downstream  recirculation  region.  Dennis  and  Smith  (1980) 
admit  that  they  do  not  resolve  this  feature  of  the  flow  and  that  the  use  of  very 
fine  grids  is  necessary  to  recover  the  true  behaviour.  However,  from  their 
numerical  calculations  they  are  able  to  imply  that  the  fluid  separates  downstream  of 
the  constriction. 

The  behaviour  of  the  vorticity  along  the  downstream  channel  wall  (y  =  x 
0)  indicates  whether  the  domain  has  been  truncated  far  enough  downstream.  As 
x  —  oo  the  vorticity  f  =  -v2#  ->  12.  If  f(x,  =;)  settles  down  to  this  value  then  the 

-  i 

truncation  length  is  adequate.  The  value  of  ?(x,  k)  is  plotted  in  Figs.  3(a)  —  (c) 
for  Re  =  0,  100,  500,  respectively.  These  figures  suggest  that  the  exit  length  is 
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suitable  for  the  range  of  Reynolds  numbers  considered  here.  Dennis  and  Smith 
(1980)  use  downstream  truncation  lengths  of  2  and  2.5  and  obtain  excellent  agreement 
in  their  results  which  demonstrates  that  the  numerical  solutions  obtained  with  a 
truncation  length  of  2  are  satisfactory.  Hunt  (1989)  uses  a  transformation  of  the 
independent  variables  to  set  the  downstream  boundary  at  x  ~  1000.  Instead  of 
experimenting  with  truncation  lengths  he  needs  to  determine  appropriate  mapping 
parameters  which  appear  in  the  transformation. 

Contours  of  the  stream  function  are  presented  in  Fig.  4(a)  —  (0  for  Re  =  0, 
10,  50,  100,  150,  200,  300,  400  and  500,  respectively,  in  the  region  —  1  <;  x  £  1, 
0  ^  y  <;  1.  The  vortex  in  the  salient  corner  diminishes  in  size  from  Re-0  until 
around  Re  =  45  and  then  increases  slowly  as  Re  is  increased  further.  Let  Lv 
denote  the  distance  between  the  point  where  the  separation  line  meets  the  top  of 
the  channel  and  the  salient  corner.  The  value  of  Lv  is  recorded  in  Table  I  for 
different  methods  and  for  a  range  of  values  of  Re.  The  results  of  Dennis  and 
Smith  (1980)  are  obtained  by  means  of  two  successive  h2  —  extrapolation  operations 
on  information  calculated  on  grids  of  mesh  lengths  1/10,  1/20  and  1/40.  Hunt  (1989) 
uses  a  transformed  grid  with  48  x  128  points.  The  values  of  Lv  in  columns  (b)  and 
(c)  agree  to  within  5%.  For  the  range  of  values  of  Re  considered  here  the  scheme 
of  Hunt  (1989)  that  uses  artificial  viscosity  gives  results  that  are  closer  to  those  in 
column  (b)  than  the  scheme  without  artificial  viscosity. 

Close-up  views,  —5  ^  x  <;  0,  5  <£  y  ^  1,  of  the  salient  corner  are  given  in 
o  o 

Fig.  5(a)  —  (e)  for  Re  =  10,  50,  100,  150  and  200,  respectively.  The  first  closed 
streamline  corresponds  to  a  contour  height  of  1  +  10~5  and  the  remaining  ones  differ 
by  multiples  of  10  s.  The  interesting  feature  in  these  plots  is  the  second  vortex 
appearing  close  to  the  corner.  In  the  Stokes  case  Moffatt  (1964)  predicts  an  infinite 
sequence  of  eddies  running  into  the  corner.  The  size  of  this  second  vortex  grows 
moderately  as  Re  is  increased.  Dennis  and  Smith  (1980)  only  just  detect  the  second 
vortex  at  values  of  Re  of  at  least  a  thousand  and  then  only  on  extremely  fine 
meshes  of  size  1/80. 

A  small  recirculation  region  downstream  of  the  constriction  is  observed  in 
Fig.  4  (d)  for  Re  =  100.  This  region  suddenly  grows  when  Re  is  increased  to  200. 
In  fact  when  the  stream  function  is  calculated  at  intermediate  values  of  Re  this 
downstream  recirculation  region  grows  suddenly  at  a  value  of  Re  around  175.  A 
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magnification  of  this  region  is  shown  in  Fig.  6  (a)  -  (d)  for  Re  -  200,  300,  400  and 
500,  respectively.  In  these  figures  the  stream  function  is  contoured  in  the 
rectangle  0  <;  x  <,  1,  |  <£  y  <;  i.  These  figures  give  a  good  description  of  the  way 
in  which  the  vortex  develops  as  Re  is  increased.  In  Fig.  6  the  fluid  is  seen  to 
separate  at  the  corner  and  not  slightly  to  the  right  as  Hunt  (1989)  observes. 
Further  the  length  of  the  downstream  recirculation  regions  are  significantly  longer 
than  those  calculated  by  Hunt  (1989)  as  can  be  seen  in  Table  II.  Dennis  and  Smith 
(1980)  do  not  detect  this  region  although  they  say  that  its  existence  is  implied  as 
the  grid  is  refined.  For  Re  -  500  they  predict  that  separation  occurs  at  a  point  in 
the  interval  0  <  x  <  0.3  with  reattachment  at  a  point  beyond  x  =  1.2.  For  this 
value  of  Re  their  prediction  of  the  reattachment  point  is  closer  to  our  calculated 
value  than  the  value  obtained  by  Hunt  (1989). 

7.  Concluding  Remarks 

The  steady  planar  two-dimensional  laminar  flow  of  an  incompressible  fluid 
through  an  abruptly  contracting  channel  is  considered  for  moderate  values  of  the 
Reynolds  number.  The  governing  equation  for  the  stream  function  is  linearized 
using  Newton’s  method  and  solved  using  spectral  domain  decomposition  techniques. 
A  conforming  domain  decomposition  is  used  which  in  conjunction  with  a  carefully 
constructed  collocation  strategy  ensures  that  the  resulting  spectral  approximations 
are  globally  -  continuous.  In  addition,  the  spectral  approximations  are  C°° 
except  along  subdomain  interfaces.  An  efficient  direct  method  for  solving  the 
spectral  collocation  equations  at  each  stage  of  the  Newton  process  is  described. 

At  most  six  Newton  iterations  are  required  for  convergence.  For  0  <  Re  <; 
150  a  converged  numerical  solution  is  obtained  from  an  initial  approximation  that  is 
zero  everywhere.  Continuation  in  Re  in  increments  of  50  is  used  for 
150  <  Re  <;500.  The  vortex  in  the  salient  corner  decreases  in  size  from  Re  -  0  to 
around  Re  -  45  and  then  grows  slowly  as  Re  is  increased  further.  A  small 
recirculation  region  just  downstream  of  the  constriction  appears  at  Re  -  100.  This 
region  suddenly  grows  as  Re  is  increased  to  a  value  around  175.  This  recirculation 
region  extends  further  downstream  as  Re  is  increased. 

There  is  qualitative  and  quantitative  agreement  with  Dennis  and  Smith  (1980) 
in  the  bulk  of  the  flow  and  the  description  of  the  salient  corner  vortex.  Their 
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numerical  solutions  do  not  predict  a  downstream  recirculation  region  although  they 
say  there  is  a  hint  of  its  existence  as  the  grid  is  refined.  Our  calculations  predict 
that  the  fluid  separates  at  the  reentrant  corner  not  just  to  the  right  of  it  as  does 
Hunt  (1989)  and  that  reattachment  occurs  further  downstream  than  in  his 
simulations. 

This  paper  demonstrates  that  spectral  domain  decomposition  techniques  are 
capable  of  solving  complex  flow  situations  and  resolving  the  mam  features  of  the 
flow.  This  is  accomplished  in  an  efficient  manner  using  far  fewer  degrees  of 
freedom  than  other  methods  of  discretization. 
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TABLE  I 


Re 

(a) 

(b) 

(c) 

(d) 

(e) 

0 

0.285 

0.290 

0.284 

- 

- 

10 

9.150 

0.148 

0.155 

- 

- 

25 

0.129 

0.125 

- 

- 

- 

50 

0.129 

0.123 

0.129 

- 

- 

75 

0.135 

0.135 

- 

- 

- 

100 

0.143 

0.140 

0.144 

- 

- 

125 

0.154 

0.150 

- 

0.164 

0.168 

150 

0.160 

0.155 

- 

- 

- 

200 

0.183 

0.183 

- 

- 

- 

250 

0.210 

0.205 

- 

0.209 

0.227 

300 

0.223 

0.223 

- 

- 

- 

350 

0.235 

0.233 

- 

- 

- 

400 

0.244 

0.244 

- 

- 

- 

450 

0.251 

0.25b 

- 

- 

- 

500 

0.260 

0.265 

0.266 

0.260 

0.308 

The  length 

of  the  upstream 

vortex  as  a 

function  of 

Re  for  (a)  spectral 

collocation 

method  on  grid  (3),  (b)  spectral  collocation  method  on  grid  (4),  (c)  extrapolated  finite 
difference  scheme  of  Dennis  and  Smith  (1980),  (d)  finite  difference  scheme  of  Hunt 
(1989)  with  artificial  viscosity,  (e)  finite  difference  scheme  of  Hunt  (1989)  without 
artificial  viscosity. 
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TABLE  II 


Re 

(a) 

(b) 

(c) 

150 

0.088 

- 

- 

:oo 

0.363 

- 

- 

250 

0.500 

0.106 

0.096 

300 

0.538 

- 

- 

350 

0.700 

- 

- 

400 

0.738 

- 

- 

450 

0.924 

- 

- 

500 

0.995 

0.406 

0.406 

The  length  of  the  downstream  vortex  as  a  function  of  Re  for  (a)  spectral 
collocation  method  on  grid  (4),  (b)  finite  difference  scheme  of  Hunt  (1989)  with 
artificial  viscosity,  (c)  finite  difference  scheme  of  Hunt  (1989)  without  artificial 
viscosity. 


CAPTIONS  FOR  FIGURES 

Fig.  3.  Wall  vortcity  f(x,l/2)  for  x  >  0,  for  various  Reynolds  numbers  Re. 

Fig.  4.  Streamf unction  contours  in  -1  ^  x  1,  0  <;  y  ^  1,  for  various  Reynolds 
numbers  Re. 

Fig.  5.  Close-up  views  of  the  salient  corner  streamlines  on  -1/8  <;  x  <,  0, 

7/8  ^  y  ^  1  for  various  Reynolds  numbers  Re. 

i 

Fig.  6.  Close-up  views  of  the  streamlines  in  the  region  0  x  <;  1,  £  ^  y  ^  ^ 
for  various  Reynolds  numbers  Re. 
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